#define _SIZE_T_DEFINED
#ifndef __CUDACC__
#define __CUDACC__
#endif
#ifndef __cplusplus
#define __cplusplus
#endif

#include <cuda.h>
#include <device_launch_parameters.h>
#include <texture_fetch_functions.h>
#include "float.h"
#include <builtin_types.h>
#include <vector_functions.h>

static __inline__ __host__ __device__ float4 fminf(float4 a, float4 b)
{
	return make_float4(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z), fminf(a.w,b.w));
}
static __inline__ __host__ __device__ float4 fmaxf(float4 a, float4 b)
{
	return make_float4(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z), fmaxf(a.w,b.w));
}

extern "C"  
{
#define FILTER_TABLE_SIZE 16
#define Gaussian2x2_xWidth 2.f
#define Gaussian2x2_yWidth 2.f
#define Gaussian2x2_invXWidth (1.f / Gaussian2x2_xWidth)
#define Gaussian2x2_invYWidth (1.f / Gaussian2x2_invXWidth)
#define GAMMA_TABLE_SIZE 1024u
#define PARAM_IMAGE_WIDTH 1280
#define	PARAM_IMAGE_HEIGHT 720
#define PARAM_IMAGE_FILTER_WIDTH_X 1.2f
#define PARAM_IMAGE_FILTER_WIDTH_Y 1.2f
#define PARAM_IMAGE_FILTER_MITCHELL_B 1 / 3.f;
#define PARAM_IMAGE_FILTER_MITCHELL_C 1 / 3.f;

	typedef unsigned int uint;

	__device__ void PixelIndex2XY(const uint index, uint *x, uint *y) {
		*y = index / PARAM_IMAGE_WIDTH;
		*x = index - (*y) * PARAM_IMAGE_WIDTH;
	}

	__device__ uint XY2PixelIndex(const uint x, const uint y) {
		return x + y * PARAM_IMAGE_WIDTH;
	}

	__device__ uint XY2FrameBufferIndex(const int x, const int y) {
		return x + 1 + (y + 1) * (PARAM_IMAGE_WIDTH + 2);
	}

	__device__ bool IsValidPixelXY(const int x, const int y) {
		return (x >= 0) && (x < PARAM_IMAGE_WIDTH) && (y >= 0) && (y < PARAM_IMAGE_HEIGHT);
	} 


	__device__ float Mitchell1D(float x) {
		const float B = PARAM_IMAGE_FILTER_MITCHELL_B;
		const float C = PARAM_IMAGE_FILTER_MITCHELL_C;

		if (x >= 1.f)
			return 0.f;
		x = fabs(2.f * x);

		if (x > 1.f)
			return (((-B / 6.f - C) * x + (B + 5.f * C)) * x +
			(-2.f * B - 8.f * C)) * x + (4.f / 3.f * B + 4.f * C);
		else
			return ((2.f - 1.5f * B - C) * x +
			(-3.f + 2.f * B + C)) * x * x +
			(1.f - B / 3.f);
	}

	// Mitchell Filter
	__device__ float ImageFilter_Evaluate(const float x, const float y) {
		const float distance = sqrt(
			x * x * (1.f / (PARAM_IMAGE_FILTER_WIDTH_X * PARAM_IMAGE_FILTER_WIDTH_X)) +
			y * y * (1.f / (PARAM_IMAGE_FILTER_WIDTH_Y * PARAM_IMAGE_FILTER_WIDTH_Y)));

		return Mitchell1D(distance);
	}

	typedef struct {
		float c1, c2, c3;
	} Spectrum;

	typedef Spectrum Pixel; 

	typedef struct {
		Spectrum radiance;
		float weight;
	} SamplePixel;

	typedef struct {
		float screenX, screenY;
		float PathWeight;
		int pass;
		Spectrum radiance;
	} SampleBufferElem;

	int Ceil2Int(const float val) {
		return (int)ceil(val);
	}

	int Floor2Int(const float val) {
		return (int)floor(val);	
	}


	__device__ float Clamp(float val, float low, float high) {
		return (val > low) ? ((val < high) ? val : high) : low;
	}

	__device__ unsigned int Floor2UInt(const float val) {
		return (val > 0.f) ? ((unsigned int)floor(val)) : 0;
	}

	__device__ float Radiance2PixelFloat(
		const float x,
		const float *gammaTable) {
			//return powf(Clamp(x, 0.f, 1.f), 1.f / 2.2f);

			const unsigned int index = min(
				Floor2UInt(GAMMA_TABLE_SIZE * Clamp(x, 0.f, 1.f)),
				GAMMA_TABLE_SIZE - 1u);
			return gammaTable[index];
	}

	__device__ float4 cross(const float4 &a, const float4 &b) {
		return make_float4
			(a.y * b.z - a.z * b.y,
			a.z * b.x - a.x * b.z,
			a.x * b.y - a.y * b.x,
			0.f);
	}

	__device__ float dot(const float4 &a, const float4 &b)
	{
		return a.x*b.x+ a.y*b.y+a.z*b.z +a.w*b.w;
	}

	__device__ float4 operator+(const float4 &a, const float4 &b) {

		return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
	}
	__device__ float4 operator*(const float4 &a, const float4 &b) {

		return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
	}
	__device__ float4 operator+=(const float4 &a, const float4 &b) {
		return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
	}
	__device__ float4 operator/(const float &a, const float4 &b) {

		return make_float4(a / b.x, a / b.y, a / b.z, 0.f);
	}
	__device__ float4 operator-(const float4 &a, const float4 &b) {

		return make_float4(a.x-b.x,a.y-b.y,a.z-b.z,a.w-b.w);
	}

	__device__ void AddWeightedSample(SamplePixel *sp, const float4 sample) {
		float4 weight = make_float4(sample.w, sample.w, sample.w, 1.f);
		float4 *p = ( float4 *)sp;

		*p += weight * sample;
	} 

	__device__ inline void AtomicAdd(float* address, float value)
	{
		float old = value;  
		float new_old;

		do
		{
			new_old = atomicExch(address, 0.0f);
			new_old += old;
		}
		while ((old = atomicExch(address, new_old))!=0.0f);

	}

	inline __device__ void AddSample1( SamplePixel *sp, const float4 sample) {

		AtomicAdd(&sp->radiance.c1, sample.x);
		AtomicAdd(&sp->radiance.c2, sample.y);
		AtomicAdd(&sp->radiance.c3, sample.z);
		AtomicAdd(&sp->weight,      sample.w);
	} 




	__device__ void AddSample( SamplePixel *sp, const float4 sample) {
		float4 *p = (float4 *)sp;
		*p += sample;
	} 

	//work_group_size(64, 1, 1)
	__global__ void PixelAddSampleBuffer(
		const unsigned int width,
		const unsigned int height,
		SamplePixel *sampleFrameBuffer,
		const unsigned int sampleCount,
		SampleBufferElem *sampleBuff) 
	{
		const unsigned int index = blockIdx.x*blockDim.x+threadIdx.x;
		if (index >= sampleCount)
			return;

		SampleBufferElem *sampleElem = &sampleBuff[index];
		const unsigned int x = (unsigned int)sampleElem->screenX;
		const unsigned int y = (unsigned int)sampleElem->screenY;
		const float4 sample = make_float4(sampleElem->radiance.c1, sampleElem->radiance.c2, sampleElem->radiance.c3, 1.f);

		AddSample1(&sampleFrameBuffer[x + y * width], sample);
		/*
		atomicAdd(&sampleFrameBuffer[x + y * width].radiance.c1, sampleElem->radiance.c1);
		atomicAdd(&sampleFrameBuffer[x + y * width].radiance.c2, sampleElem->radiance.c2);
		atomicAdd(&sampleFrameBuffer[x + y * width].radiance.c3, sampleElem->radiance.c3);
		atomicAdd(&sampleFrameBuffer[x + y * width].weight, sampleElem->PathWeight);
		*/
		/*

		sampleFrameBuffer[x + y * width].radiance.c1 += sampleElem->radiance.c1;
		sampleFrameBuffer[x + y * width].radiance.c2 += sampleElem->radiance.c2;
		sampleFrameBuffer[x + y * width].radiance.c3 += sampleElem->radiance.c3;
		sampleFrameBuffer[x + y * width].weight +=1.f;
		*/
		//AddWeightedSample(&sampleFrameBuffer[x + y * width], sample);
		//AddSample(&sampleFrameBuffer[x + y * width], sample);
	}

	//work group size 8, 8, 1
	__global__ void PixelClearFB(
		const unsigned int width,
		const unsigned int height,
		Pixel *frameBuffer) 
	{
		const unsigned int px = blockIdx.x*blockDim.x+threadIdx.x;;
		if(px >= width)
			return;
		const unsigned int py = blockIdx.y*blockDim.y+threadIdx.y;
		if(py >= height)
			return;
		const unsigned int offset = px + py * width;

		Pixel *p = &frameBuffer[offset];
		p->c1 = 0.f;
		p->c2 = 0.f;
		p->c3 = 0.f;
	} 
	//work group size 8, 8, 1

	__global__ void PixelClearSampleFB(
		const unsigned int width,
		const unsigned int height,
		SamplePixel *sampleFrameBuffer) {
			const unsigned int px = blockIdx.x*blockDim.x+threadIdx.x;;
			if(px >= width)
				return;
			const unsigned int py = blockIdx.y*blockDim.y+threadIdx.y;
			if(py >= height)
				return;
			const unsigned int offset = px + py * width;

			SamplePixel *sp = &sampleFrameBuffer[offset];
			sp->radiance.c1 = 0.f;
			sp->radiance.c2 = 0.f;
			sp->radiance.c3 = 0.f;
			sp->weight = 0.f;
	}

	//work group size 8, 8, 1

	__global__ void PixelUpdateFrameBuffer(
		const unsigned int width,
		const unsigned int height,
		SamplePixel *sampleFrameBuffer,
		Pixel* frameBuffer,
		const float* gammaTable) {
			const unsigned int px = blockIdx.x*blockDim.x+threadIdx.x;;
			if(px >= width)
				return;
			const unsigned int py = blockIdx.y*blockDim.y+threadIdx.y;
			if(py >= height)
				return;
			const unsigned int offset = px + py * width;

			SamplePixel *sp = &sampleFrameBuffer[offset];
			Pixel *p = &frameBuffer[offset];

			const float weight = sp->weight;
			if (weight == 0.f)
				return;

			const float invWeight = 1.9f / weight;
			p->c1 = Radiance2PixelFloat(sp->radiance.c1 * invWeight, gammaTable);
			p->c2 = Radiance2PixelFloat(sp->radiance.c2 * invWeight, gammaTable);			
			p->c3 = Radiance2PixelFloat(sp->radiance.c3 * invWeight, gammaTable);
	} 
}